Graph Theory and Kidney Transplants
This post will introduce some concepts from the mathematical discipline of graph theory, and go on to explain how some of the techniques developed in this field can be used to help facilitate kidney exchanges. I was inspired to write about this topic by the chapter ‘The Kidneys of Königsberg’ in Ian Stewart’s popular mathematics book What’s the Use?, and I thoroughly recommend both that chapter, and the rest of the book!
Some Graph Theory Background
Suppose that you are a teacher in charge of taking a group of six students, who we will call Students A, B, C, D, E, and F, on a school trip, and you need to decide what pairs they are going to sit in on the minibus. You ask the students to anonymously let you know who they would be willing to sit with and find out that:
- Student A and Student D would be happy to be in a pair.
- Student A and Student E would be happy to be in a pair.
- Student A and Student F would be happy to be in a pair.
- Student B and Student C would be happy to be in a pair.
- Student B and Student E would be happy to be in a pair.
- Student C and Student D would be happy to be in a pair.
We can represent this information in a graph, which, in discrete mathematics, means a structure consisting of vertices and edges (as opposed to describing a chart or a plot).
Here, the points of the graph – their vertices (sometimes also called nodes) – represent the students, and the lines between the vertices, known as edges, encode the possible pairings they can sit in; if two students would accept being seated together, they have an edge between them. Graphs are commonly used to study community structures in sociology, and they can also be applied to protein-protein interaction networks in biology and transportation networks in urban planning, to give just two other examples.
Graphs can be represented computationally in several different ways; here, we will discuss three of the most common – adjacency matrices, adjacency lists, and edge lists. For a graph \(G\) with \(n\) vertices, the adjacency matrix of \(G\) is an \(n\times n\) matrix \(M\), where the entry in row \(i\) and column \(j\) of \(M\) is \(1\) if vertex \(i\) and vertex \(j\) are connected by an edge, and \(0\) if they are not. For the rest of this post, we will denote the specific graph structure above, representing the possible seating arrangements of the students, by \(S\). We can define the adjacency matrix for \(\text{S}\) in R as follows:
S_adj_mat <- matrix(c(# Let's walk through how we fill out the first row of the adjacency matrix:
# Row 1 represents the people Student A can sit with;
# they form a compatible pair with Students D, E, and F,
# so there is a 1 in column 4, representing Student D,
# a 1 in column 5, representing Student E,
# and a 1 in column 6, representing Student F.
# The pairings A-B and A-C are unacceptable so there are 0s in columns 2 and 3.
# Vertex A is not connected to itself, because Student A cannot form a pair with themselves,
# so there is a 0 in column 1 – in fact, all the diagonal values of this adjacency matrix are 0
# because none of the vertices have edges connecting them to themselves.
0, 0, 0, 1, 1, 1,
0, 0, 1, 0, 1, 0,
0, 1, 0, 1, 0, 0,
1, 0, 1, 0, 0, 0,
1, 1, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0),
nrow = 6, byrow = TRUE)
print(S_adj_mat)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 1 1 1
## [2,] 0 0 1 0 1 0
## [3,] 0 1 0 1 0 0
## [4,] 1 0 1 0 0 0
## [5,] 1 1 0 0 0 0
## [6,] 1 0 0 0 0 0As well as an adjacency matrix \(M\) for \(G\), we can also create an adjacency list, \(L_A\), which is a list of length \(n\), where each entry is also a list. The \(i\)th element of \(L_A\) is a list of vertices that vertex \(i\) of graph \(G\) forms an edge with – this list can have a length of \(0\) if vertex \(i\) is not connected to any vertices.
# Let's write a function which takes as its input a graph's adjacency matrix, and returns the adjacency list for the graph
adj_mat_to_adj_list <- function(adj_mat){
num_nodes <- ncol(adj_mat) # we find the number of vertices of the graph
adjacency_list <- vector(mode = "list", length = num_nodes) # we initialise our adjacency list as an empty list with the same number of elements of vertices of the graph
# Now we make an empty list for each vertex and fill it with the indexes of the vertices that it is connected to:
for (i in 1:num_nodes){
connected_vertices <- list()
for (j in 1:num_nodes) {
if (adj_mat[i,j] == 1) {
connected_vertices <- c(connected_vertices, j)
}
}
adjacency_list[[i]] <- connected_vertices
}
adjacency_list
}
# We can use this function to find the adjacency list for S:
S_adj_list <- adj_mat_to_adj_list(S_adj_mat)
names(S_adj_list) <- c("A","B","C","D","E","F")
# And find all the vertices that are connected to vertex A, for example:
vert_A_connections <- S_adj_list$A
LETTERS[unlist(vert_A_connections)]## [1] "D" "E" "F"Finally, we can also create an ‘edge list’, \(L_E\), of all the edges of \(G\); for the graph we have been looking at
this would be A-D, A-E, A-F, B-C, B-E, C-D. In R, we can
represent each edge as a vector of length two where the elements are the
vertices that the edge is connected to, and thus encode \(L_E\) as a two column matrix.
# We can write a simple function to get the edge list matrix of a graph from it's adjacency matrix:
adj_mat_to_edge_list <- function(adj_mat) {
edge_list <- unique(t(apply(which(adj_mat == 1, arr.ind = TRUE), 1, sort)))
edge_list
}
# Let's see what the output looks like for S:
S_edge_list <- adj_mat_to_edge_list(S_adj_mat)
print(adj_mat_to_edge_list(S_adj_mat))## [,1] [,2]
## [1,] 1 4
## [2,] 1 5
## [3,] 1 6
## [4,] 2 3
## [5,] 2 5
## [6,] 3 4# We can convert the vertex numbers to letters to check we've got the result we've expected:
apply(S_edge_list, 1:2, function(i) LETTERS[i])## [,1] [,2]
## [1,] "A" "D"
## [2,] "A" "E"
## [3,] "A" "F"
## [4,] "B" "C"
## [5,] "B" "E"
## [6,] "C" "D"Now that we know some definitions, we can return to our problem –
finding a configuration in which the maximum possible number of
passengers on the minibus are happy with their partner. Each person can
only sit with one other person, so in graph theory terms we are looking
to find a matching of the graph – this is a set of
edges of the graph in which no two edges share any common vertices1. But it’s
not enough just to find any matching – for example the set of edges
\(S_{1}\)=A-D is a
matching, as is the set \(S_{2}=\)C-D, but neither are
very useful to us as they both only give us one pairing.
Instead, we are looking for a maximum matching, which is a matching of a graph that has the largest possible number of edges. One way to find a maximum matching of a graph \(G\) would be to check all the possible subsets of edges of \(G\), throw out the ones that aren’t matchings, and then pick the largest remaining subset. But since the set of all possible subsets of edges is equal to the power set of the set of edges of \(G\), it contains \(2^{E}\) subsets for us to check, where \(E\) is the number of edges of \(G\). For our graph \(S\), this would give us \(2^6=64\) subsets to check, which might just about be manageable, but as the number of edges increases, the number of subsets we will need to check using this naive approach will grow exponentially. A graph with 12 edges, for example, would require us to check \(2^{12}=4096\) subsets!
xs <- 1:12
ys <- 2^xs
plot(xs, ys, pch=4, cex=1.5, lwd=3, col="#009ed8",
xlab="Number of Edges of Graph", ylab="Number of Subsets to Check")
To understand the algorithms that allow us to find a maximum matching
more quickly, we need another concept from graph theory –
augmenting paths. Given a graph \(G\), and a matching \(M\), an augmenting path is a sequence of
edges, \(A\), that connect two
unmatched vertices of \(G\) such that
the edges \(A\) are alternately in
\(M\) and not in \(M\). For example, if we have graph with
four vertices and an edge list 1-2, 2-3, 3-4 and we are
given a matching2-3, the augmenting path is
1-2, 2-3, 3-4.
Below, I’ve written (not very well/efficiently!) an R function that
will find an augmenting path for a given graph and matching. It takes as
input the edge list of the graph, and a vector, match_vec
which represents the matching. The length of match_vec
should be the same as the number of vertices of the graph, and
match_vec[i] is the index of the vertex that the \(i\)th vertex is matched to – if it is
unmatched then we’ll define match_vec[i] to be zero.
augmenting_path <- function(adj_list, match_vec) {
vertices_match_status <- match_vec != 0 # we create a logical vector to represent whether each vertex is included in the matching
num_vertices <- length(vertices_match_status) # we find the number of vertices in the graph
unmatched_vertices <- which(vertices_match_status == FALSE, arr.ind=TRUE) # we create a vector of the vertices which are not in the matching
for (i in 1:num_vertices) { # we will try finding an augmenting path starting from each unmatched vertex until we find one
aug_path <- matrix(, nrow=0, ncol=2) # we initialise an empty matrix to start the augmenting path
if (vertices_match_status[i] == TRUE) { # we skip the matched vertices
next
}
root <- i
current_vert <- i # we will use this to keep track of which vertex we have reached in our attempt to find an augmenting path
# Let's keep track of the vertices we've included in our augmenting path attempt from this vertex:
unvisited_from_root <- c(1:num_vertices)
unvisited_from_root <- unvisited_from_root[!unvisited_from_root==root]
repeat { # we create a repeating loop that will run the same procedure until we break out of it using break or return
neighbours <- unlist(adj_list[[current_vert]]) # we create a vector of the vertices that have an edge between them and the vertex that we are currently on
neighbours_unvisited_from_root <- intersect(neighbours, unvisited_from_root) # we create a vector of vertices adjacent to the vertex we are currently on that have not already been included in our attempt to find an augmenting path
if (length(neighbours_unvisited_from_root) == 0) { # if there are none of these, we can stop running our repear loop and try our searh again from the next unmatched vertex
break
}
match_to_current_vert <- match_vec[current_vert]
# if we are on a vertex that is part of the matching, and the vertex it is matched to is not yet part of the augmenting path, we want to go to that vertex next
if (match_to_current_vert %in% unvisited_from_root) {
next_vert <- match_to_current_vert
unvisited_from_root <- unvisited_from_root[!unvisited_from_root==next_vert]
aug_path <- rbind(aug_path, c(current_vert,next_vert))
current_vert <- next_vert
next
}
unmatched_neighbours_unvisited_from_root <- intersect(unmatched_vertices, neighbours_unvisited_from_root) # we create a vector of unmatched vertices adjacent to the vertex we are currently on that we haven't already included in our augmenting path attempt
if (length(unmatched_neighbours_unvisited_from_root) > 0) { # if any of these exist, then we have found an augmenting path, and we can return it!
aug_path <- rbind(aug_path, c(current_vert,unmatched_neighbours_unvisited_from_root[1]))
return(aug_path)
}
# If we didn't return an augmenting path in the if statement above, all the vertices adjacent to the vertex we are currently on that are not already included in our augmenting path attempt must be part of the matching, so we will go to one of these next and continue searching:
next_vert <- neighbours_unvisited_from_root[1]
unvisited_from_root <- unvisited_from_root[!unvisited_from_root==next_vert]
aug_path <- rbind(aug_path, c(current_vert,next_vert))
current_vert <- next_vert
}
}
# If we haven't returned an augmenting path after searching from all unmatched vertices, we return an empty matrix:
aug_path <- matrix(, nrow=0, ncol=2)
return(aug_path)
}Let’s check to see if this function works on our toy example of a graph with four vertices and a one edge matching.
example_adj_mat <- matrix(c(0,1,0,0,
1,0,1,0,
0,1,0,1,
0,0,1,0), nrow = 4, byrow = TRUE)
example_adj_list <- adj_mat_to_adj_list(example_adj_mat)
example_match <- c(0,3,2,0)
augmenting_path(example_adj_list, example_match)## [,1] [,2]
## [1,] 1 2
## [2,] 2 3
## [3,] 3 4So my augmenting_path function is not very pretty, but
it does the trick! Now, how will this help us find the maximum matching
of a graph? Well in 1957, mathematician Claude Berge proved that if an
augmenting path cannot be found for a graph and a matching, then that
matching is a maximum matching2. This is the insight that is central to the
Hungarian
Maximum Matching Algorithm, which works by starting with any
matching on a graph (including an empty matching) and searching the
graph from each unmatched vertex until it finds an augmenting path. It
then finds the symmetric
difference of the path and the current matching, and this set of
edges becomes the new, improved matching that it will use for its next
iteration. For example, in the scenario above where we have a graph with
edge list 1-2, 2-3, 3-4 and a matching 2-3,
the first iteration of the algorithm will find the augmenting path
1-2, 2-3, 3-4, and will update the matching to be
1-2, 3-4.
When it can no longer find any augmenting paths, the Hungarian Maximum Matching Algorithm terminates, and thanks to Claude Berge, we can be sure that the matching it has found is a maximum matching.
# Let's use our augmenting_path function to write a Hungarian Maximum Matching Algorithm that takes as its input the adjacency matrix of a graph and a matching vector:
maximum_matching <- function(adj_mat, match_vec) {
adj_list <- adj_mat_to_adj_list(adj_mat) # we use the function we wrote earlian
num_vertices <- length(match_vec) # we find the number of vertices in the graph
# Let's create a matrix to store the current matching:
current_matching <- matrix(, nrow=0, ncol=2)
# We can populate it with the initial matching to start with:
for (i in 1:num_vertices) {
if (match_vec[i] == 0) {
next
}
current_matching <- rbind(current_matching, c(i,match_vec[i]))
}
if (nrow(current_matching) != 0) {
current_matching <- unique(t(apply(current_matching, 1, sort))) # we order the edges with the lower indexed vertex first, and remove duplicates
}
# Now we will run our augmenting path algorithm and and use it to increase our current matching until no more augmenting paths can be found:
current_match_vec <- match_vec
repeat {
aug_path <- augmenting_path(adj_list, current_match_vec) # we find an augmenting path for the input graph and the current matching
aug_path_length <- nrow(aug_path)
# If the augmenting_path function has not found any augmenting paths, we break out of the repeat loop:
if (aug_path_length == 0) {
break
}
aug_path <- t(apply(aug_path, 1, sort)) # we order the edges in the augmenting path with the lowest indexed vertex first
# If the augmenting path consists of only one edge, we add that edge to our matching and run another iteration of the augmenting_path function with the new matching:
if (aug_path_length == 1) {
current_matching <- rbind(current_matching, aug_path)
# Before running the augmenting_path function again we need to update the current_match_vec algorithm to reflect our new matching:
current_matching_length <- nrow(current_matching)
current_match_vec <- rep(0,num_vertices)
for (i in 1:current_matching_length) {
current_match_vec[current_matching[i,1]] <- current_matching[i,2]
current_match_vec[current_matching[i,2]] <- current_matching[i,1]
}
next
}
# If the augmenting path has more than one edge, we add the edges that are not included in the current matching, and remove the ones that are
# (by the definition of an augmenting path, these alternate):
for (i in 1:aug_path_length) {
if (i%%2 == 1) {
current_matching <- rbind(current_matching, aug_path[i,])
}
if (i%%2 == 0) {
current_matching <- current_matching[!(current_matching[,1]==aug_path[i,1] & current_matching[,2]==aug_path[i,2]),]
}
}
# Before running the augmenting_path function again we need to update the current_match_vec algorithm to reflect our new matching:
current_matching_length <- nrow(current_matching)
current_match_vec <- rep(0,num_vertices)
for (i in 1:current_matching_length) {
current_match_vec[current_matching[i,1]] <- current_matching[i,2]
current_match_vec[current_matching[i,2]] <- current_matching[i,1]
}
}
return(current_matching[order(current_matching[,1]),]) # we return the first matching we arrived at for which no augmenting paths could be found,
# with the edges ordered by their lowest indexed vertex
}What does this algorithm give us when applied to our example graph \(S\), and an empty matching?
S_empty_matching <- rep(0,6)
S_max_matching <- maximum_matching(S_adj_mat, S_empty_matching)
apply(S_max_matching, 1:2, function(i) LETTERS[i])## [,1] [,2]
## [1,] "A" "F"
## [2,] "B" "E"
## [3,] "C" "D"So students A and F, B and E, and C and D, can all sit in pairs.
This may seem like something we could have worked out easily by hand, but matching algorithms are really useful when there are so many vertices and edges to a graph that a human could never work out a maximum matching just from looking at it.
You may be relieved to learn that we are nearly at the end of the background graph theory I wanted to include in this article before discussing matching algorithms’ applications to kidney transplants! However, before we move on, I should admit that I have been oversimplifying things somewhat. Though the Hungarian Maximum Matching Algorithm will work on a huge variety of graphs, there are some graphs that will trip it up. For example, consider the graph and matching below, taken from this stack exchange post by Misha Lavrov.
Trying to find an augmenting path from vertex 1, the algorithm will
go from 1 to 2 to 8 to 9 to 10, before getting stuck and concluding that
there is no augmenting path, and thus that the given matching is a
maximum matching. Similarly, if it starts from vertex 4 it will follow
the sequence of edges 4-3,3-7,7-6,6-5 before it again gets
stuck and terminates.
Graphs that run into these kinds of problems are always non-bipartite, meaning they contain cycles of odd numbers of edges, like the pentagons in the graph above. To ensure that, if they exist, augmenting paths can always be found in these types of graph, the Hungarian Maximum Matching Algorithm needs to be modified. The first person to find a way to do this was computer scientist Jack Edmonds, who introduced an ingenious method for adapting the algorithm to deal with odd length cycles in his ‘blossom’ algorithm in 19653. The details of his algorithm are fascinating, but beyond the scope of this post, which is why I have only used bipartite graphs as examples, and will continue to do so in the next section. If you want to learn how a maximum matching can be found for any graph using Edmonds’ blossom algorithm, I recommend this excellent youtube video by Tomáš Sláma.
What Does Graph Theory Have to Do With Kidney Donations?
People with end-stage renal disease have kidneys that are unable to perform their function of filtering the blood and are, therefore, unable to survive without undergoing regular dialysis unless they receive a working kidney transplant from a donor. Many people with this condition will have family member or loved one who is willing to donate a kidney to them, as most of us only need one working kidney to live a healthy life, but this is often not possible due to a mismatch between the donor and the intended recipient’s blood and/or antigen type. For example, imagine that a patient with advanced chronic kidney disease, who we’ll call Recipient A, has a sibling, Donor A, who is willing to donate a kidney to them but cannot because they have incompatible blood groups. In the past, this would have meant that Donor A would have been unable to donate a kidney in order to help their sibling4. Nowadays, however, the pair have another option; they can participate in a the UK Living Kidney Sharing Scheme5. Suppose another recipient and donor pair, a parent and child, say, are in the same situation. If the donor in this pair – let’s call them Donor F for reasons that will shortly become clear – is a match with Recipient A, and Donor A’s kidney is compatible with the recipient, Recipient F, then the two pairs can effectively ‘swap’ kidneys, in what is called a paired exchange.
We can represent possible paired kidney exchanges between pairs of
donors and recipients in a graph, with each vertex representing one
incompatible pair, and an edge linking two pairs signifying that a
paired exchange is possible between them. Consider the scenario
described above; this would be represented by a graph with two vertices
and one edge, A-B.
What happens if we introduce more incompatible pairs of donors and recipients? Well, we might end up with a graph like \(S\), but this time with the vertices representing patients in need of kidneys and loved ones who are unable to denote to them. In reality a graph of incompatible donor-recipient pairs in a country like the UK would have hundreds of vertices, but we can use a smaller one to help us think about the problem of matching compatible pairs.
As we saw earlier, several different matchings of this graph exist, but some of them will result in more people receiving kidneys than others. For example, if the Pair A’s transplant centre search a regional or national database for a pair that Pair can exchange a kidney with and stop at the first one they find, then Pair A might be matched with Pair D. This would leave Pair B and Pair C free to participate in a paired exchange, but Recipients E and F would then be left without a kidney. If a maximum matching algorithm is used on the other hand, all six recipients will be matched with a compatible donor.
In 2004, a team from John Hopkins Hospital in the US proposed using Edmonds’ blossom algorithm to optimize kidney exchanges6. They used simulated pools of donor-recipient pairs to test their scheme, and found that it would result in more transplants, with better antigen matches, than the ‘first-accept’ system that was in use at the time7. Of course, employing an algorithm to decide who gets which organs presents numerous ethical and logistical challenges. For instance, what happens when there is more than one possible maximum matching, resulting in different people getting kidney transplants? Consider what our example graph would look like if we discovered that Recipient B was actually not compatible with Donor E.
